function Fisher = mvnrfish(Data, Design, Covar, MatrixFormat, CovarFormat)
%MVNRFISH Fisher information for multivariate normal or least-squares regression.
%	Fisher information matrix based on current maximum likelihood or least-squares parameter
%	estimates.
%
%		Fisher = mvnrfish(Data, Design, Covar);
%		Fisher = mvnrfish(Data, Design, Covar, MatrixFormat);
%
% Inputs:
%	Data - NUMSAMPLES x NUMSERIES matrix with NUMSAMPLES samples of a NUMSERIES-dimensional random
%		vector. If a data sample has missing values (with NaNs), the sample is ignored.
%	Design - Either a matrix or a cell-array to handle two distinct model structures. First, if
%		NUMSERIES = 1, Design can be a NUMSAMPLES x NUMPARAMS matrix with known values. This is the
%		"standard" form for regression on a single data series. Alternatively, for any
%		NUMSERIES >= 1, Design can be a cell array of either 1 or NUMSAMPLES cells, where each cell
%		contains a NUMSERIES x NUMPARAMS matrix of known values. If Design has a single cell, then
%		it is assumed to be the same Design matrix for each sample. Otherwise, Design must contain
%		individual Design matrices for each sample.
%	Covar - NUMSERIES x NUMSERIES matrix of estimates for the covariance of the residuals of the
%		regression.
%
% Optional Inputs:
%	MatrixFormat - String that identifies parameters to be included in the Fisher information
%		matrix. The default method is 'full'. The choices are:
%		'full' - (default) Compute the full Fisher information matrix.
%		'paramonly' - Compute only components of the Fisher information matrix associated with the
%			model parameter estimates.
%	CovarFormat - String that specifies the format for the input covariance matrix. The choices are:
%		'full' - Default method. Compute the full covariance matrix.
%		'diagonal' - Treat the covariance matrix as a diagonal matrix.
%
% Outputs:
%   Fisher - TOTALPARAMS x TOTALPARAMS Fisher information matrix based on current parameter
%		estimates, where
%			TOTALPARAMS = NUMPARAMS + NUMSERIES * (NUMSERIES + 1)/2
%		if MatrixFormat = 'full' and
%			TOTALPARAMS = NUMPARAMS
%		if MatrixFormat = 'paramonly'.
%
% WARNING: If calculating the full Fisher information matrix, this function is VERY slow.
%
% See also MVNRSTD, MVNRMLE.

%	Copyright 2005-2007 The MathWorks, Inc.
%	$Revision: 1.1.6.2 $ $Date: 2007/05/10 13:45:02 $

% Step 1 - check arguments

if nargin < 5 || isempty(CovarFormat)
	CovarFormat = 'FULL';
end
if nargin < 4 || isempty(MatrixFormat)
	MatrixFormat = 'FULL';
end
if nargin < 3
	error('Finance:mvnrfish:MissingInputArg', ...
		'Missing required input arguments Data, Design, or Covar.');
end
if isempty(Data)
	error('Finance:mvnrfish:EmptyDataArray', ...
		'The required input argument Data is empty.');
end
if isempty(Design)
	error('Finance:mvnrfish:EmptyDesignArray', ...
		'The required input argument Design is empty.');
end
if isempty(Covar)
	error('Finance:mvnrfish:EmptyCovar', ...
		'The required input argument Covar is empty.');
end

%[NumSamples, NumSeries, NumParams] = ...
%	checkmvnrsetup(Data, Design, [], Covar, true);

[NumSamples, NumSeries] = size(Data);
if iscell(Design)
	if (numel(Design) == 1)
		SingleDesign = true;
	else
		SingleDesign = false;
	end
	NumParams = size(Design{1},2);
else
	SingleDesign = false;
	NumParams = size(Design,2);
end

if ~all(size(Covar) == [NumSeries, NumSeries])
	error('Finance:mvnrfish:InconsistentDims', ...
		'The covariance matrix Covar has wrong dimensions.');
else
	[CholCovar, CholState] = chol(Covar);
	if CholState > 0
		error('Finance:mvnrfish:NonPosDefCov', ...
			'Covariance matrix is not positive-definite.');
	end
end

% Step 2 - initialization

if ~any(strcmpi(MatrixFormat,{'PARAMONLY','FULL'}))
	warning('Finance:mvnrfish:UnknownFormatString', ...
		'Unknown MatrixFormat string. Will use default FULL.');
	MatrixFormat = 'FULL';
end

if strcmpi(MatrixFormat,'PARAMONLY')
	TotalParams = NumParams;
else
	if strcmpi(CovarFormat,'FULL')
		TotalParams = NumParams + (NumSeries * (NumSeries + 1))/2;
	else
		TotalParams = NumParams + NumSeries;
	end
end

Fisher = zeros(TotalParams,TotalParams);

% Step 3 - do partials wrt Mean

Count = 0;
TestMatrix = zeros(NumParams,NumParams);

if iscell(Design)
	if SingleDesign
		A = CholCovar' \ Design{1};
		for k = 1:NumSamples
			TestMatrix = TestMatrix + A'*A;
		end
		Count = NumSamples;
	else
		for k = 1:NumSamples
			if ~any(isnan(Data(k,:)))
				Count = Count + 1;
				A = CholCovar' \ Design{k};
				TestMatrix = TestMatrix + A'*A;
			end
		end
	end
else
	for k = 1:NumSamples
		if ~any(isnan(Data(k,:)))
			Count = Count + 1;
			A = CholCovar' \ Design(k,:);
			TestMatrix = TestMatrix + A'*A;
		end
	end
end
TestMatrix = (1.0/Count) .* TestMatrix;

for i = 1:NumParams
    for j = 1:i
        Fisher(i,j) = TestMatrix(i,j);
        Fisher(j,i) = Fisher(i,j);
    end
end

% Step 4 - do partials wrt Covar

if strcmpi(MatrixFormat,'FULL')
	if strcmpi(CovarFormat,'FULL')
		InvCovar = inv(Covar);

		GradC1 = zeros(NumSeries,NumSeries);
		GradC2 = zeros(NumSeries,NumSeries);

		i = NumParams;
		for i1 = 1:NumSeries
			for j1 = 1:i1
				i = i + 1;

				GradC1(i1,j1) = 1;						% do dC/dtheta(i)
				GradC1(j1,i1) = 1;

				j = NumParams;
				for i2 = 1:NumSeries
					for j2 = 1:i2
						j = j + 1;

						if (j <= i)
							GradC2(i2,j2) = 1;			% do dC/dtheta(j)
							GradC2(j2,i2) = 1;

							Temp1 = InvCovar*GradC1;
							Temp2 = InvCovar*GradC2;

							Fisher(i,j) = 0.5*trace(Temp1*Temp2);
							Fisher(j,i) = Fisher(i,j);                        

							GradC2(i2,j2) = 0;			% undo dC/dtheta(j)
							GradC2(j2,i2) = 0;
						end
					end
				end

				GradC1(i1,j1) = 0;						% undo dC/dtheta(i)
				GradC1(j1,i1) = 0;
			end
		end
	else		% diagonal covariance
		CDiag = 1 ./ diag(Covar);
		FDiag = 0.5 * CDiag .^2;
		Fisher(NumParams*TotalParams + NumParams+1:TotalParams+1:end) = FDiag;
	end
end
